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A recent experiment has considered the effective permeability of two-phase flow of air and a 
water-glycerol solution under steady-state conditions in a two-dimensional model porous medium, 
and found a power law dependence with respect to capillary number. Running simulations on a 
two-dimensional network model similar power law behavior is found, for high viscosity contrast as 
in the experiment and also for viscosity matched fluids. Two states are found, one with stagnant 
clusters and one without. For the stagnant cluster state, a power law exponent 0.50 is found for 
viscosity matched fluids and 0.54 for large viscosity contrast. 

PACS numbers: 47.56.+r, 89.75.Da, 47.61. Jd, 47.55.Ca 



I. INTRODUCTION 

A system consisting of two immiscible fluids flowing 
simultaneously through a porous medium has provided a 
complex challenge for both scientists and engineers [l|, [|J ■ 
An improved understanding of two-phase flow in porous 
media is of considerable importance for such applications 
as enhanced oil recovery and groundwater management. 
Also, attention has recently been directed towards the re- 
duced efficiency of polymer electrolyte membrane (PEM) 
fuel cells caused by produced water flooding the pore 
space of the PEM's gas diffusion and catalyst layers ||. 

Simulations of two-phase flow in porous media may be 
done by modeling the porous medium as a network |^-|8| 
or by simulating fluid dynamics directly on porous me- 
dia reconstructions from microstructure images [9]. The 
advantage of the latter approach is a more realistic de- 
scription of the pore space. Network modeling, on the 
other hand, has the advantage of allowing a large range of 
simulation parameters on reasonably large systems, even 
when the viscous pressure field needs to be re-solved for 
each timestep. 

Most work on two-phase flow in porous media, both nu- 
merical and experimental, has focused on displacement 
processes [10|. Relative permeabilities, used in reservoir 
modeling, can be measured under such un-steady con- 
ditions. Another approach is to maintain constant ex- 
ternal control parameters, such as total flow-rate and 
average saturation, and run the simulation until macro- 
scopic quantities such as relative permeabilities only fluc- 
tuate around mean values. These are steady-state values, 
which in general will be different from the un-steady val- 
ues. An open question is whether there exists some range 
of control parameters within which the steady-state val- 
ues are unique, independent of initial conditions. 

Network steady-state simulations have considered vari- 
ous aspects such as the dynamics of disconnected ganglia 
of oil, the phase diagram in the regime where both phases 
are mobile, the transition where one of the phases be- 
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comes immobilized and diffusion on non- wetting clusters 

[mil. 

Experimental work on steady-state with simultaneous 
flow of both phases has been done in two dimensions us- 
ing chamber-throat networks etched in glass [l9l - l23| and 
Hele-Shaw cells [24], HH . A recent experiment performed 
with three-dimensional bead packings [2(| is also of in- 
terest to this work. Recently, microfluidic devices have 
been used as model porous media, achieving unprece- 
dented control over the pore-scale geometry, but so far 
they have only been used to study drainage [271 ] . 

When studying steady-state two parameters are usu- 
ally kept under control. Flow is then maintained until all 
measurable quantities only fluctuate around some mean 
value. One control parameter is the capillary number Ca, 
a dimensionless quantity which parametrizes the relative 
strength of viscous and capillary pressure drops. It is 
proportional to the total flow-rate. The other parameter 
which controls the steady-state may be either the non- 
wetting saturation S or the fractional flow of the non- 
wetting phase F. The pressure drop required to drive 
flow may be measured, as well as relative permeabili- 
ties. The sum of relative permeabilites is the effective 
permeability K c g. This is typically found to be less than 
unity under steady-state conditions, which signifies that 
the mixture of two phases results in a lower permeability 
than if only a single phase is present. 

Tallakstad et al. [H, [25[ have reported that the pres- 
sure drop AP driving flow of air and a water-glycerol so- 
lution under steady-state conditions in a two-dimensional 
model porous medium shows a power law dependence 
with respect to the capillary number Ca, 

AP-Ca' 3 . (1) 

They keep fractional flow fixed at F = 7/15, due to in- 
jecting the two phases through 15 syringes of equal flow- 
rates, 7 of which inject air. In the following this will be 
referred to as the Oslo experiment. 

Rassi et al. (2(| have also reported a power law similar 
to Eq. [TJ Their experiment uses three-dimensional bead 
packings partially saturated with air under flowing wa- 
ter conditions. The partially saturated state is achieved 
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by initial simultaneous flow of both two phases, combin- 
ing them at a tee junction beneath the bead pack. In 
the following this experiment will be referred to as the 
Bozeman experiment. 

Using a two-dimensional network simulation this work 
investigates whether power law behavior similar to Eq. [1] 
is found when saturation is kept fixed, and whether it ex- 
tends beyond the limit of large viscosity ratios used in the 
experiments. A preliminary report has been published in 



II. THE MODEL 
A. General 

The porous medium is modeled as a network of con- 
nected capillaries. A capillary corresponds to the throat 
connecting two pores, within which one or more menisci 
may be present. With regards to the pressure drop across 
a meniscus the throats are modeled as having an hour 
glass shape, while the flow equations are solved with per- 
meabilities corresponding to cylindrical tubes. This is 
a matter of convenience, and reflects an interest in the 
network-scale dynamics rather than a realistic descrip- 
tion at the pore scale. 

The network considered in this work is a square lat- 
tice inclined at 45° relative to the main direction of flow. 
Boundary conditions are bi-periodic, such that the net- 
work may be mapped onto a torus. Disorder is intro- 
duced by assigning random radii r to throats, from a flat 
distribution between 10% and 40% of the throat length. 
All throats are assumed to be equally long, with length 
I. 

A modified version of the Young-Laplace relation is 
used to obtain the capillary pressure drop 



2a f 2ttx. 
Pc(x) = — I 1 - cos (— . 



(2) 



where x is the position of the meniscus in the interval 
[0,£] and a is the surface tension. 

Disregarding wetting films and considering flow to be 
laminar, the Hagen-Poiseuille permeability for cylindrical 
tubes gives the flow-rate q by the Washburn equation, 
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where \i is the viscosity of the phase contained in the 
throat and Ap is the difference between the pressures of 
the two nodes connected by the throat. If the viscosities 
of the two phases are different, the viscosity /z of Eq. [3] is 
replaced by the volume-weighted effective viscosity aver- 
aged over the length of the throat. The two phases are 
considered to be incompressible. 

Flow is driven by a pressure drop applied across a cut 
through the network. The cut is invisible to the local 
flow field, so it does not cause any boundary effects. A 



constant volumetric flow-rate Q is achieved by solving 
the pressure field for two applied pressure drops AP and 
extrapolating linearly to the pressure drop which gives 
the desired Q using the relationship 



Q = Qo + < 



kq <f>A 

Meff L y 



AP + Q c , 



(4) 



where the first term is the single-phase flow-rate Qo cor- 
responding to Darcy's law and Q c is a flow-rate induced 
by the capillary pressure drops, i.e., it is the flow-rate 
the system would have if the external pressure drop was 
set to zero. Under steady-state conditions Q c is always 
in the opposite direction of Q, i.e., the net effect of the 
capillary pressure drops is to oppose flow. In the Darcy 
term Qo, L y is the length of the sample in the direction of 
flow, A is the cross-sectional area of the sample and <f> is 
its porosity, kq is the single-phase absolute permeability 
of the porous medium and fi e R is the volume-weighted 
effective viscosity averaged over the entire system. 

After solving the pressure field, the simulation is iter- 
ated forward in time. The throat with largest flow-rate is 
identified and the timestep is set by allowing any menisci 
in this throat to advance 10% of the throat length. Most 
throats will have flow-rates orders of magnitude slower 
than this, giving a smooth motion of menisci through 
the hour glass shaped throats. The 10% rule is chosen by 
trial and error as being as large as possible without affect- 
ing the results of the simulation. Timestepping is done 
by the Euler method, which is only first order accurate. 
The numerical error may be estimated by comparing the 
difference between driving power and heat dissipation, 
which should average out to zero in the steady-state. It 
is found that the error is negligible for Ca > 10~ 3 , but 
grows to a few percent for Ca = 10~ 4 . For Ca lower than 
this the Euler method is therefore not reliable. 

Since the network is a closed system (mapped onto a 
torus) the total volume of the two phases is conserved, 
i.e., any given run of the simulation maintains a constant 
saturation S which is determined by the initialization 
conditions. 

When a meniscus reaches the end of a throat, new 
menisci are created in throats flowing out of the con- 
necting pore. In order to prevent the number of menisci 
from growing without bonds some scheme is needed to 
merge nearby pairs of menisci. Menisci within a certain 
distance from the ends of throats may be merged, mim- 
icking fluid coalescence withing the pore volumes. The 
closest pair of menisci may be merged whenever some 
maximum allowed number of menisci within a throat is 
reached. A minimum separation between menisci may 
also be enforced, merging pairs of menisci when they 
come too close. Whenever pairs of menisci are merged, 
this is done in such a way that the merging does not 
contribute to fractional flow. All these schemes require 
some coalescence parameter. All have been implemented 
and tested. It is found that while fractional flow depends 
on the merging scheme and coalescence parameter, the 
power law exponent does not. Only the range of validity 
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of the power law is affected, as it ceases to be valid when 
fractional flow becomes either or f . 

For the results given here, pairs of menisci less than 
10% of the length of a throat apart are merged. 

A more realistic description would need to carefully 
consider mechanisms for snap-off and coalescence. As a 
detailed description of flow at the pore scale is not the 
primary interest of this work, the simpler approach is 
considered sufficient. 

An impression of the model is provided by Fig. [1] Con- 
nected strands of large non-wetting saturations can be 
seen (dark) amidst patches of high wetting saturations 
(light). Internal flow-rates span several orders of mag- 
nitude. Channels of both high and low flow-rates are 
clearly visible. There are channels within channels, form- 
ing a complex pattern. For the choice of parameters and 
initialization shown in Fig. [1] channels are dynamic. At 
any given time flow-rates outside fast channels are rel- 
atively slow, but when menisci block old channels the 
blocks are by-passed by new channels, and no spatial re- 
gion remains stagnant indefinitely. 



B. Definitions 

The main parameter of the model is the relative 
strength of capillary and viscous forces, parametrized by 
the capillary number 



Ca: 



<j(j)A ' 



(5) 



Larger values of Ca correspond to faster flow, or lower 
surface tension. 

The viscosity ratio M is defined as 



M = 



(0) 



where subscripts denote the non-wetting and wetting 
phases. In this work a viscosity ratio M = 1 (viscos- 
ity matched fluids) and M = 10~ 4 (gas/liquid) as in the 
experiment has been explored. 

The effective permeability re e ff is defined by 



Keff 



which gives 



Q_ 
Qo' 



k cS k 4>A 
Mcff L y 



AP. 



(7) 



(8) 



Eq. [8] resembles Darcy's law, with K e sKo substituted for 
the absolute permeability. It should be kept in mind 
that K e g depends on Q and S, so Eq. [8] should not be 
interpreted as a linear relation between Q and AP. 

A power law dependence between effective permeabil- 
ity and capillary number is expressed as 



which is meant to be valid within some range for Ca. 

Given Eq. |9] and considering that Q ~ Ca when all 
other parameters are fixed, the scaling of AP can be 
expressed as 



AP-K^Q-Ca 1 " 7 , 



(10) 



K ff ~ Ca 7 , 



(9) 



which together with the definition in Eq. [T] gives the rela- 
tion 7 + /3 = 1. Eqs.fTlandfTUlare valid when Q ~ Ca, i.e., 
all parameters in Ca except Q are kept constant. This is 
the case in the Oslo and Bozeman experiments. Eq. pre- 
lates two dimensionless quantities and is valid regardless 
of the manner by which Ca is changed. 



C. Initialization and steady-state 

The reached state is not necessarily independent of how 
the initial conditions are chosen - different initial distri- 
butions of the two phases may lead to different states. 
For this reason it is important to compare different ini- 
tialization schemes. 

In order to achieve a desired saturation the throats 
are initially filled with only one phase, after which the 
other phase is inserted in randomly chosen throats until 
the desired saturation is reached. Two schemes for in- 
serting a phase into a throat are considered here. First, 
the chosen throat may be filled completely. Second, a 
bubble of some small, fixed volume and consisting of two 
menisci may be inserted at some random position within 
the throat. Several bubbles may then exist in a single 
throat. If they overlap, they are merged. In the first 
scheme there are no capillary pressure drops initially, as 
these vanish at the ends of the throats. In the second 
scheme the total capillary pressure drop across the sys- 
tem averages out to zero, but there are finite, random 
capillary pressure drops on the scale of the pore throats. 
These two schemes will be referred to as large bubble and 
small bubble initialization, respectively. 

After the desired saturation is reached, two schemes 
are considered for the approach to steady-state. First, 
the simulation may be run with a surface tension which 
is initially low, and then gradually increased until the 
desired Ca is reached. Second, the simulation may be 
run with constant parameters throughout. This will be 
referred to as gradual and abrupt initialization, respec- 
tively. 

After a sufficient number of timesteps the system al- 
ways settles into a dynamic steady-state, unless a span- 
ning cluster of one of the phases is formed, in which case 
the spanning cluster carries all flow and the system ceases 
to fluctuate. This will be referred to as the spanning clus- 
ter transition. The dynamic steady-state is characterized 
by persistent fluctuations around average values with no 
perceivable drift. It is observed that the range of param- 
eters which lead to a dynamic steady-state coincides with 
those which gives flow of both phases. 
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FIG. 1. Snapshot of the spatial distribution of local saturations s (left) and absolute-valued flow-rates \q\ (right) under 
steady-state conditions. The flow-rates are plotted on a logarithmic scale, in arbitrary units. The overall direction of flow is 
up. Axes are measured in lattice units (£). Bi-periodic boundary conditions are used. M = 1, Ca = 10 -3 and S = 0.5. 



III. RESULTS 

2 m 1 1 1 r 



A. Two states 

It is found that there are two possible states, one state 
with stagnant clusters and one without. 

The stagnant cluster state is reached by large bubble 
and abrupt initialization, for a narrow range of satura- 
tions S ~ 0.2. The stagnant clusters are spatial regions 
with slow flow-rates and slowly changing saturation dis- 
tributions. The locations of the stagnant clusters are 
static throughout the simulations. The stagnant clus- 
ters form after only a very small number of timesteps 
and their locations depend on the initial distribution of 
phases. In this sense, the stagnant cluster state is a state 
with memory of initial conditions. The meaning of this 
is discussed in Sec. IIVB1 The range of saturations which 
lead to the stagnant cluster state is affected by the choice 
of merging scheme and coalescence parameter. 

Reaching the stagnant cluster state requires particular 
choices for initialization scheme and parameters. By con- 
trast, the state without stagnant clusters is reached by all 
initialization schemes and parameters, except those that 
produce the stagnant cluster state or a spanning cluster. 
No memory of initial conditions is observed. There are 
structures: channels of flow and connected strands of the 
non-wetting phase, but these do not appear to have any 
particular length-scale. All throats and menisci appear 
to affect and be affected by the overall flow. This is in 
contrast with the stagnant cluster state, where the config- 
urations within stagnant clusters appear to be insulated 
from the overall flow pattern. Because all menisci par- 
ticipate in forming the pattern of flow, the state without 
stagnant clusters will be referred to as the foam state. 

Starting from the stagnant cluster state, the stag- 
nant clusters will mobilize if Ca is increased above some 




FIG. 2. (Color online) Scaling of K e ff with Ca for saturations 
0.1 (circles), 0.2 (triangles) and 0.4 (squares). Straight lines 
are power law fits below a threshold and equal to 1 above 
it. Error bars are omitted as they are smaller than the point 
representations. M = 1. System size is 32x64. 



threshold. If Ca is then returned to its original value, 
the foam state results. The stagnant cluster state can 
thus transition to the foam state given a perturbation. 
The opposite transition is never observed. This leads to 
identifying the stagnant cluster state as metastable. The 
foam state, on the other hand, is robust with regards to 
changes in Ca - after a period of increased or decreased 
Ca the system returns to the same state. 
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B. Power law exponents 

Fig. [2] shows simulation results for K e g in the foam 
state, for a broad range of Ca and M = 1. Three val- 
ues of S are shown. Four simulation runs underlie each 
data point. Each run has independent realizations of the 
network disorder and initial phase distribution. The sim- 
ulation is ended when the initial transient is judged to 
occupy no more than the initial one-third of the run; K e g 
is then averaged over the last third. The variation of the 
averaged K e ff between independent runs is quite small. If 
the variations were plotted as error bars on Fig. [21 they 
would be smaller than the point representations. 

A power law is fitted by eye, with a transition at some 
upper threshold value for Ca, around which the system 
has a smooth but distinct transition. The threshold value 
is denoted Ca up . The fit for n c g is given by 



Keff = < 



Ca 
Ca, 



up 



Ca <c Ca U p, 
Ca ^> Ca u p. 



(11) 



From Fig. [5] a distinct bending is observed with respect 
to the fit. A possible cause of the bending may be the nu- 
merical error caused by the Euler timestepping method, 
which becomes more serious at low Ca. The bending may 
also signal some cross-over behavior between two regimes 
with different scaling exponents. Also, if the power law 
is governed by some critical point at a finite Ca, omit- 
ting this will cause an apparent bending near the critical 
point. A further discussion of a possible modification of 
the power law is given in Sec. IIVCI 

Fig. [3] shows simulation results for K g in the stagnant 
cluster state. Some bending exists in the limit of very 
small Ca. It is in this limit that error bars begin to 
become substantial. For a range of values of Ca exceeding 
two decades, good power law behavior is observed. 

Tab. U reports the best fit values for 7 and Ca up for 
the foam and stagnant cluster states, for M = I and 
M = 10~ 4 and various values of S. In the foam state the 
best fit for the exponent 7 depends on S 1 , and varies in 
the range of 0.67 - 0.80. 

At the large Ca limit of the law (Ca > Ca up ) the per- 
meability approaches the single-phase permeability of the 
system, K e g — 1. The existence of this viscous transition 
is not surprising; it indicates the point where capillary 
pressure drops become perturbations of a dominating 
field of viscous pressure drops. The transition has been 
investigated for hysteresis, which has not been found. Its 
smooth, reversible behavior is reminiscent of a second 
order phase transition. For Ca > Ca up , F = S for all 
values of S. A precise estimation of Ca up has not been 
attempted. The results presented in Tab. Q] are those re- 
sulting from the fit according to Eq. [TT] 

The power law also has a limit for low Ca. When Ca is 
decreased below some threshold Cai , a spanning cluster 
forms. At this spanning cluster transition the saturation 
distribution and flow pattern becomes static. Decreasing 




FIG. 3. (Color online) Scaling of n e g with Ca for S = 0.2 in 
the stagnant cluster state. Below the threshold Ca up , a power 
law with 7 = 0.50 is shown. Error bars are included, as these 
become larger than the point representations for Ca < 10~ 4 . 
M = 1. System size is 32x64. 



TABLE I. The scaling exponent 7 and threshold value Ca up 
for different values of M and S. System size is 32x64. 



M 



7 



Ca u 



Stagnant cluster state: 



1 

10" 4 
Foam state: 

1 

1 

1 

1 
10 
10 
10 
10 



0.2 
0.2 

0.1 
0.2 
0.3 
0.4 
0.1 
0.2 
0.3 
0.4 



0.50 
0.54 

0.68 
0.72 
0.76 
0.80 
0.67 
0.71 
0.75 
0.78 



0.042 
0.040 

0.015 
0.030 
0.035 
0.040 
0.017 
0.029 
0.037 
0.040 



Ca further after the formation of a spanning cluster does 
not change K c g. Independent runs at Ca beneath the 
spanning cluster transition display significant hysteresis, 
and the concept of a reproducible state is lost. Due to 
the hysteresis effects, the spanning cluster transition is 
reminiscent of a first order phase transition. A further 
discussion of this is given in [l7| . 

Caio depends on S. For the purpose of this work, the 
transition is avoided by choosing S in the range of 0.1 — 
0.4. For larger or smaller values of S the spanning cluster 
transition would appear within the range of Ca shown in 
Fig. El 

No size-dependence has been found for K c g. In addi- 
tion to the results for 32x64 systems in Tab. |H identical 
values have been found for 16x32 systems. Also, sev- 
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eral representative points have been checked for larger 
systems, up to 256x512, with indistinguishable results. 

IV. DISCUSSION 
A. The viscosity ratio 

It is somewhat surprising that K e & has almost no de- 
pendence on the viscosity ratio M. The simulations dis- 
play a strong dependence on M whenever transients are 
considered, in agreement with established experimental 
results. To the author's knowledge, no experiment has 
studied a dependence on M under steady-state condi- 
tions. 

The reason for the lack of sensitivity to M shown by 
the simulations may be attributed to the lack of large 
single-phase clusters under steady-state conditions. If 
the distribution of menisci causes the system to resemble 
a foam, no flow of one phase is possible without causing 
the other phase to flow in nearly equal measure. Then, 
only the volume-weighted effective viscosity fj, c g matters; 
the viscosity ratio M does not. 

In the high Ca limit this explanation seems plausible. 
It is in this limit that the picture of a foam with small 
variations in the saturation distribution is most appro- 
priate. However, when Ca is reduced the distribution 
of saturation becomes more heterogeneous, and a large 
viscosity contrast could be expected to become more im- 
portant. The absence of an increased M dependence at 
low Ca indicates that the relative importance of capillary 
pressure drops relative to viscous pressure drops, which 
increases with decreasing Ca, out-weighs the increased 
importance of M. 

In other words, the scaling of K e g with Ca is dominated 
by the effects of capillary pressure drops, with only a 
modest contribution from heterogeneities in the viscosity 
distribution. 

The viscosity ratio has a significant importance in at 
least one respect. The spanning cluster transition at Cai 
is sensitive to M. For M>1 (more-viscous non-wetting 
fluid), Caio is between 10 -2 and 10~ 3 for all values of S. 
This does not allow for a power law spanning much more 
than a single decade. For this reason the case of M> 1 
is not considered in this work. 



B. The metastable state 

How can the stagnant cluster state be claimed to be 
a state, given that it only appears for a particular ini- 
tialization scheme and the location of stagnant clusters 
depends on the initial distribution of phases? While the 
state retains memory of the initial distribution of phases 
for a long time, possibly indefinitely, the macroscopic av- 
erages are consistent from realization to realization. The 
fluctuations of effective permeability within a single run 
and the fluctuations of run averages between independent 



runs are of the same order. This is true both when the 
individual runs are on identical networks with different 
initial phase distributions and when the individual runs 
have different realizations of the network disorder. The 
variation of macroscopic averages between different sim- 
ulation runs are all consistent with these fluctuations. In 
this sense, the macroscopic averages are not sensitive to 
the initial phase distribution. 

In other words, the stagnant cluster state has well- 
defined properties despite displaying hysteresis. For this 
reason, it is referred to as a state. 

It should be noted that a stagnant cluster does not 
consist of a single phase. Due to the way the simulations 
are initialized, both phases will in general exist within the 
cluster. Therefore, there will be some capillary pressure 
drops within. These pressure drops do not seem to be or- 
ganized, they are random and cancel out. This indicates 
that the perimeter of the stagnant cluster is formed early 
in the simulation, before the menisci within the cluster 
organize themselves in response to the flow. 

If S is larger than 0.3 the stagnant clusters are always 
mobilized, eventually resulting in the foam state, even for 
the case of large bubble and abrupt initialization. This is 
due to the larger pressure drop required to obtain the de- 
sired flow-rate at these saturations. The larger pressure 
drop prevents stagnant clusters from forming, and the 
foam state results. For a discussion on the dependence 
of AP on S, see 

C. The scaling of the foam state 

A recent paper (30j explores a power law which is dif- 
ferent from Eq. |H1 The simulator used in their paper is 
similar to the one used here. They explore a power law 
on the form (AP - AP C ) ~ Ca^ 3 and find /3 w 0.5. AP C is 
a threshold pressure drop, beneath which the authors of 
[30| state that no flow occurs. The introduction of this 
new power law is in part motivated by the wandering 
value of (3 reported in [28[ . 

The primary purpose of this work is to consider the 
power law given by Eq. |H1 It corresponds to the power 
laws reported by the Oslo and Bozeman experiments. 
However, the bending observed in Fig. [5] does indicate 
that a power law of a different form should also be con- 
sidered. 

In order to explore this possibility, a more general 
power law is considered in the following, on the form 

(k cS - K e ff,cr) ~ (Ca - Ca cr ) 71 , (12) 

where Cai Q < Ca < Ca up . Ca cr and K e ff,cr = K cff 
may now be interpreted as critical values for the control 
and order parameter, respectively. In the following they 
will be treated as independent parameters in an opti- 
mization problem, aiming to obtain the best possible fit 
to the power law. 

Part of the motivation for considering Eq. Q2] is to see 
whether a dependence of Ca cr on S can remove the appar- 



7 




io -4 io -3 ^o 2 io _1 10° 
Ca 



FIG. 4. (Color online) Scaling of K e ff with Ca for saturations 
0.2 (triangles) and 0.4 (squares). The porous medium was 
initialized without disorder. Otherwise the simulations are 
identical to those in Fig. H M = 1. S = 0.2 has 7 = 0.82 and 
Ca up = 0.025. S = 0.4 has 7 = 0.85 and Ca up = 0.048. 
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FIG. 5. (Color online) Scaling of k b $ — K e ff,cr with Ca — Ca cr 
for saturations 0.1 (circles), 0.2 (triangles) and 0.4 (squares). 
Straight lines are power law fits. M = 1. 



ent dependence of the exponent on S. If an ^-dependent 
critical value is ignored this can be misinterpreted as an 
S'-dependent exponent [3l| . 

In addition to the simulation results reported for the 
foam state in Tab. HI simulations without disorder are 
considered here. Fig. @] shows results from simulations 
with M = 1, S = 0.2 and S = 0.4, identical to those 
reported in Fig. [3] except that all throat radii are set 
equal to r — 0.21. With no disorder, the exponent 7 
increases to 0.82 and 0.85 and less bending is observed. 
As before, the power laws have been fitted by eye. 

The critical values Ca cr and K e s,cr are determined by 
minimizing the error of the best fit for the power law. 
Instead of fitting by eye, the best fit is now found using 
a least squares method, considering only the data points 
within the interval 10~ 4 < Ca < 10~ 2 . The error is 
minimized with repsect to variations of the two critical 
values by means of simulated annealing. 

Tab. HTlreports the results from this procedure. Fig. [5] 
displays the results and power law fits for the same values 
of M and S as Fig. H 

The values found for Ca cr are beneath the range of 
Ca which are accessible to the simulator at this time. 
Therefore, the results reported in Tab. |TT] may not be 
reliable, as they are extrapolated by an assumption of a 
power law on the form of Eq. [12] and can not be verified 
independently of this. 

For the case of no disorder, the datasets give reasonable 
fits to a power law without the introduction of critical 
values. When an attempt is made to estimate the critical 
values, they are found to be orders of magnitude smaller 
than in the disorderd case. This makes extrapolating 
their values highly inaccurate, therefore they have simply 



TABLE II. The scaling exponent 71, Ca cr and « e fl,cr f° r dif- 
ferent values of M and S, with and without disorder. Only 
the foam state is considered. 



M 


S 


71 


C&cr 




rveff ,cr 


Disorder: 












1 


0.1 


0.65 


7.0 ■ 10" 


5 


1.3 ■ 10" 


1 


0.2 


0.71 


5.9 ■ 10" 


-5 


4.4 • 10' 


1 


0.3 


0.75 


6.1 ■ 10" 


5 


3.1 • 10- 


1 


0.4 


0.77 


4.0 ■ 10" 


-5 


9.8 • 10- 


10" 4 


0.1 


0.66 


7.7-10" 


5 


1.4 • 10- 


10" 4 


0.2 


0.71 


6.1 ■ 10" 


5 


5.1 • 10- 


HT 4 


0.3 


0.75 


5.2 ■ 10" 


-5 


2.4 • 10" 


10" 4 


0.4 


0.77 


3.9 ■ 10" 


5 


1.6 ■ lO- 


No disorder 












1 


0.2 


0.81 


0.0 




Cl. 


1 


0.4 


0.85 


0.0 




0.0 



been set to zero. 

The conclusion that can be drawn from this analysis, 
is that a power law on the form of Eq.[T2]can explain the 
observed bending with respect to Eq. [SJ This does not 
remove the dependence of the exponent on S. The esti- 
mated values for Ca cr are of the same order of magnitude 
as Cai . Hysteresis and numerical inaccuracies prevent 
the hypothesis Cai = Ca cr from being explored further. 
It has also been found that the exponent is sensitive to 
disorder, as are the hypothesized critical values. Fur- 
ther work should investigate whether this remains true 
in three dimensions. 
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D. Comparison with experiments — « fr(Ca) 

The Oslo experiment [H|, HH found an exponent j3 — 
0.54 ± 0.08, corresponding to j = 0.46 ± 0.07. Their 
result covers roughly two decades in Ca. Their definition 
of Ca is different from Eq. [5] by the definition used here 
their experiments are carried out at Ca roughly between 
3 • 10~ 6 and 2 • 10~ 4 . 

They use a two-dimensional Hele-Shaw cell: a mono- 
layer of glass beads sandwiched between transparent 
plates. At the inlet 15 evenly spaced syringes inject the 
two phases with equal flow-rate through each syringe. 
The two phases are injected through every other syringe, 
with 7 syringes injecting air, fixing fractional flow to 
F = 7/15. The other syringes inject a water-glycerol 
solution which will simply be referred to as water in the 
following. 

Initial conditions consist of a fully water-saturated 
medium. Initially, pockets of air grow at the injection 
points. Eventually, the pockets become detached from 
the injection points and travel through the medium until 
they escape at the far end. It is observed that smaller 
pockets are more likely to be trapped. The pockets of 
air are embedded within a spanning cluster of water. In 
the steady-state no drift is observed in the average global 
pressure drop, and the distribution of air pockets appears 
to be translationally invariant in a statistical sense. 

The Oslo experiment explores values of Ca and S 
which, interpreted in terms of the simulations, places 
them below the spanning cluster transition. In the sim- 
ulations only one of the phases flow after Ca is reduced 
below Cai . In the Oslo experiment both phases flow 
despite the existence of a spanning cluster, due to the 
simultaneous injection. Also, below the spanning cluster 
transition, hysteresis is expected. This suggests that the 
experiment only sees simultaneous flow of air and water 
due to forcing air to flow by the conditions at the inlets. 
This will in the following be referred to as forced frac- 
tional flow, which should be distinguished from the case 
of free fractional flow. Free fractional flow is the case in 
the simulations here when Ca > Cai . 

There are substantial differences between the Oslo ex- 
periment and the simulations reported here. However, 
the Oslo experiment shares some similarities with the 
stagnant cluster state. First, the stagnant clusters are 
characterized by low, nearly vanishing flow-rates sur- 
rounded by channels of fast flow, held in place by a 
perimeter of capillary pressure drops. The trapped pock- 
ets of air in the Oslo experiment must be held in place by 
a similar mechanism. Second, the stagnant clusters are 
formed by initial conditions with few capillary pressure 
drops and develop until they reach a size where they are 
by-passed by flow without being mobilized. The pockets 
of air in the Oslo experiment are formed at the inlets, 
where they are allowed to grow until they reach a size 
where flow of the surrounding water forces them to ei- 
ther mobilize or fragment. In both cases, clusters with a 
characteristic size are produced. 



Could the setup of the Oslo experiment obtain the 
foam state? Two approaches may be suggested. First, 
large values of Ca could be explored, possibly obtained by 
adding surfactants to reduce surface tension. This would 
also allow an investigation of the transition at Ca up . Sec- 
ond, instead of injecting the two phases at alternating 
syringes, the two phases could be mixed in tee junctions 
prior to injection through the syringes. This would pre- 
vent the build-up of large air pockets at the inlets, pos- 
sibly resulting in a different state. 

In [H, UK a scaling argument is presented for (3 = 
7 = 0.5. A close relative of this argument is presented in 
Sec. El where the comparison between the Oslo experi- 
ment and the stagnant cluster state is revisited. 

The Bozeman experiment [26l ] found an exponent (3 in 
the range of 0.30 — 0.45, corresponding to 7 in the range 
of 0.55 — 0.70. Their result covers a little more than a 
single decade in Ca. Their definition of Ca is similar to 
the one used here, and is roughly in the range between 
10" 4 and 10" 3 . Their definition of Ca uses the viscosity 
of water rather than the volume- weighted average of the 
two phases, however the water saturation does not change 
dramatically except for at the largest Ca. The different 
definitions do not seriously affect the exponents. 

They use a three-dimensional bead pack, where water 
floods the sample at gradually increasing flow-rates af- 
ter an initial preparation with simultaneous flow of air 
and water. The initial preparation results in a partially 
saturated state. As the sample is flooded by water, air es- 
capes until the only pockets of air left within the sample 
are immobile. Even at the largest flow-rates the pattern 
of flow is found to be different from that of the fully 
water saturated sample, meaning that there are always 
some pockets of air remaining. These pockets are sug- 
gested to be quite small, as they do not show up on MR 
images. Repeated floodings of a single sample gives com- 
parable results, suggesting that the method produces a 
reproducible state. 

The partially saturated state is achieved by a proce- 
dure which mixes the two phases before they are injected, 
as opposed to the Oslo experiment. This could prevent 
the formation of stagnant clusters. At the end of the 
initial preparation only water flows. The saturation de- 
creases gradually until the spanning cluster transition is 
reached, at which all air is trapped. When the flow-rate is 
increased, some air mobilizes temporarily until the span- 
ning cluster transition is reached again. If the experiment 
resembles the foam state as it approaches the spanning 
cluster transition, the partially saturated state will have 
an effective permeability near that of the foam state at 
similar S and Ca. 

The best estimate for /3 provided in 26] is 0.33. This 
gives 7 0.67, which could suggest that their experiment 
has indeed seen evidence of the foam state. However, 
there arc significant differences between the Bozeman 
experiment and the simulations. The simulations are 
two-dimensional and saturation is kept constant under 
steady-state conditions, as opposed to gradually decreas- 
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ing while the sample is flooded. In Sec. IV Dl it is argued 
that the stagnant cluster exponent 7 0.5 should be 
valid for both two and three-dimensional porous media. 
This is further indication the the Bozeman experiment 
does not produce the stagnant cluster state, which leaves 
the foam state. It is not known whether the foam state 
exponents change in three dimensions. 



E. Comparison with experiments — S'(Ca) 

Both the Oslo and Bozeman experiments have reported 
decreasing S with increasing Ca. In both cases it is found 
that the low Ca limit is S « 0.5. 

The Bozeman experiment finds that S strongly de- 
creases for Ca > 10 -3 . Because water is the only injected 
phase it is not surprising that S decreases for larger Ca, 
this agrees with the spanning cluster transition moving 
to lower S for larger Ca. This also suggests that the 
Bozeman experiment transitions into the regime of free 
fractional flow. If the sharp decrease is caused by a tran- 
sition from the regime of forced to free fractional flow, the 
sharp decrease can be explained as the disappearance of 
hysteresis effects. 

In the Oslo experiment the decrease of S is harder 
to explain. As Ca is increased, the characteristic size of 
trapped air pockets decreases. The number of air pockets 
must somehow be governed by the mechanisms which 
govern flow of air. More air pockets should allow more 
air to flow. The decrease of S with increasing Ca suggests 
that the number of air pockets needed to maintain a fixed 
F does not grow fast enough to compensate the reduced 
size of the air pockets. 

At the Ca up transition seen in the simulations, frac- 
tional flow curves have a distinct transition from a non- 
trivial dependence on S and Ca, to F = S for all values 
of S and all values of Ca > Ca up . This means that if the 
Oslo experiment were to increase Ca further, they should 
at some point observe a reversal in the S'(Ca) trend. S 
should increase to an ultimate value of S = F = 7/15 
above the Ca up transition. The value of Ca which gives 
a minimum value of S most likely coincides with the tran- 
sition from forced to free fractional flow. 



V. SCALING ARGUMENT FOR /3 = 7 = 1/2 

In [24], [25| an argument is presented for why /3 = 7 = 
1/2. It will be referred to as the Oslo argument in the 
following. 

An argument which is closely related to the Oslo argu- 
ment is given in Sees. IV Al and IV Bl The two arguments 
lead to the same scaling exponent. Sees. IV CI and IV Dl 
compare the two arguments and discuss their applicabil- 
ity. 



A. Postulates 

Three postulates are stated. Given that these are jus- 
tified, (9 = 7 = 1/2 results. 

The physical picture is that of two phases somehow 
distributed through and injected into a two-dimensional 
porous medium, initially with few capillary pressure 
drops. Eventually menisci with large capillary pressure 
drops develop, causing some regions to resist flow. The 
imposed flow is then carried in regions that remain open 
to flow. These two regions are the stagnant cluster vol- 
ume and the dissipative volume, respectively. The stag- 
nant clusters are held in place by a perimeter of capillary 
pressure drops and are as large as they can be without 
mobilizing. The dissipative volume consists of a con- 
nected network of channels spanning the system. 

Postulate 1. The volume V of the sample can be di- 
vided into two sub- volumes, 

V = V diss + V^tag, (13) 

where the dissipative volume Vdi ss satisfies 

diss diss 

D is the total dissipation within the sample, d is the dis- 
sipation within a throat, the sums run over the throats 
within the dissipative volume, u is the mean velocity 
within the dissipative volume and VP is the mean pres- 
sure gradient within this volume. Separating d into the 
product of the mean values of u and VP relies on flow 
within the dissipative volume respecting Darcy's law and 
being reasonably homogeneous, or fluctuations not being 
correlated. 

Postulate 2. Consider the balance of pressure gradients 
along the stagnant cluster perimeters. Everywhere along 
the perimeter there is a pressure balance p a — Pb ~ Pc, 
where p a and pb are the pressures on either side of the 
cluster perimeter and p c is the pressure drop caused by 
the meniscus on the perimeter. p c scales with the surface 
tension a, which gives 

VP- 1* oca, (15) 

where I* is the length of the path that connects the ex- 
treme points of the stagnant cluster, i.e., from the point 
of highest to the point of lowest surrounding pressure. I* 
is a characteristic length scale of the stagnant clusters, 
and should be comparable with both their mean width 
l x and mean length l y . 

Postulate 3. It is necessary to relate I* to Vdiss- As- 
sume that the mean width of the flowing channels does 
not scale with Ca. Both experiments and simulations 
observe that channels are only a single pore- width wide, 
supporting the assumption. Also, assume that the as- 
pect ratio of the stagnant clusters l y /l x does not scale 
with Ca, such that 

I* « (l X + ly) ~ l X ~ ly. (16) 
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Smaller Ca gives larger clusters, but the aspect ratio does 
not necessarily change. In experiments, clusters that mo- 
bilize are observed to be elongated along the direction of 
flow, while the smaller, stagnant clusters have l x ~ l y . 
In the simulations, the stagnant clusters appear slightly 
elongated in the direction of flow, and larger clusters 
seem to be more elongated. No Ca dependence is ob- 
served. 

The dissipative volume scales only with the number of 
channels when their width is constant. The number of 
channels A^ h scales as 



N 



eh 



L x /lx 



where L x is the width of the sample. This gives 

1 v diss' 



(17) 



(18) 



The tortuosity of the flow channels r also plays a role in 
the argument. It can be defined as the ratio of the mean 
channel length to the length of the sample L y . Given the 
assumptions above, the scaling of r with Ca is 



(l*)-H* 



const, 



(19) 



where Ny = L y /l y is the number of clusters over which 
the sum runs. Larger clusters give larger, but fewer twists 
and turns for the flow channels, which results in a tortu- 
osity r which does not scale with Ca. 



B. The scaling exponents 

Consider the effective cross-section £ e g = Vdiss /-ky 
The flow velocity must scale as 



Q Q 



Soff Vdis 



(20) 



Next, consider the scaling of the flow velocity with the 
pressure gradient, 



u ~ VP ' 



AP 



(21) 



Combining Eqs. [T5J QJ2 Wl and [2JJ gives 

Q 



u ^ VP- -~o- Vd^, (22) 

''diss ' 



which leads to 



V 2 - — ^ Ta 

K diss _ 



(23) 



Using Eg . [Til gives an expression for the total dissipation 

as 



D - VdissM ~ V di! 



Q \ 2 Q 2 Q 2 



V dl! 



Vdiss Ca v 



The dissipation can also be expressed by using Eq. [8l 

O 2 

D = APQ- — , (25) 

«eff 

from which the scaling is obtained as 

K off - Ca 1/2 =>> 7 = 1/2. (26) 

Other scaling results that can be read out of the pre- 
ceding equations are 

AP- VP- (I*)' 1 - Ca 1/2 /3 = 1/2. (27) 

In the last scaling laws, dimensional quantities are related 
to a dimensionless number. It should be understood that 
this is done for Q — Ca, which is only valid when all 
other dimensional quantities are kept constant. Eq. l26l is 
valid regardless. 

C. Comparison with the Oslo argument 

The Oslo argument is very similar to the argument 
given here. The assumptions of postulates 1, 2 and 3 
are made by both. However, there are three significant 
differences between the arguments. 

First, Eq. 16 of [25[ gives some length scale I* (not 
identically defined as the I* used here, but similar) as 
being proportional to the ratio of capillary pressure drop 
to overall viscous pressure drop. This is similar to postu- 
late 2. The argument then ignores the capillary pressure 
drop, giving I* cx AP -1 . However, the capillary pressure 
drop should not be ignored since the interest is in scal- 
ing with Ca, which contains the surface tension. Since 
the Oslo argument uses Ca — Q their conclusion is not 
affected by this, however it reduces the generality of the 
argument, which is best expressed in terms of dimension- 
less quantities. 

Second, Eqs. 16 and 17 of [25j contain the ratio AP/L, 
which is meant to be the pressure gradient surrounding 
a stagnant cluster. Taking L as the system length L y 
ignores the tortuosity of flow paths as they wrap around 
stagnant clusters. This tortuosity could in principle scale 
with Ca. However, the assumptions of postulate 3, which 
the Oslo argument also make, mean that r does not scale 
with Ca. This is an important point within the scaling 
argument and should not be ignored. If r scaled with Ca, 
the conclusion AP — Ca 1 / 2 would no longer be valid. 

Third, the Oslo argument considers the case where the 
phase contained in stagnant clusters has a very low vis- 
cosity, and presents their argument as a limiting case. 
The simulations suggest that the viscosity ratio is of lit- 
tle importance. The argument as presented here also 
does not rely on a particular viscosity ratio. If there is 
a non-unity viscosity ratio, Ca should be defined with 
the volume- weighted effective viscosity of the dissipative 
volume rather than the total volume, otherwise the argu- 
ment is unaffected by M. Any discrepancy between the 
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simulations and the argument due to different definitions 
of Ca does not affect the scaling exponents, because the 
saturation within the dissipative volume does not scale 
with Ca. The same should be true for the experiment. 



D. More on the scaling argument 

In the Oslo experiment the largest air pockets migrate 
while small ones are trapped. The migration of a large 
air pocket could lead to fragmentation. Small pockets of 
air could be released into the channels of flowing water, 
leaving behind a now trapped air pocket. The small, flow- 
ing pockets could exit the system, or coalesce with other 
trapped air pockets, causing them to migrate until they 
fragment, such that the process repeats itself. This sug- 
gests the mechanism behind forced fractional flow: over- 
sized pockets of air produced at the inlets begin to mi- 
grate, then fragment, leading to a cascade of coalescence, 
migration and fragmentation through the system. 

The migrating cluster state is a very different picture 
from that of the stagnant cluster state, where no mi- 
gration is seen. There are clearly substantial differences 
between the simulations and the Oslo experiment. How- 
ever, the origin of the /3 = 7 = 1/2 exponent appears to 
be the same in both, namely the applicability of postu- 
lates 1, 2 and 3. 

The argument is not meant to be exact. In the case 
of the Oslo experiment, the elongation of the largest air 
pockets is likely to be connected with the mechanism 
which allows fractional flow. Because the aspect ratio 
of stagnant clusters is important to postulate 3 the argu- 
ment is likely to be modified by an improved understand- 
ing of this mechanism. The separation into two volumes, 
where one contains all dissipation and obeys Darcy's law 
and the other is perfectly balanced by a perimeter of 
capillary pressure drops, is also an approximation. The 
stagnant clusters are not perfectly non-dissipative in the 
simulations, and the air pockets in the Oslo experiment 
are not all completely trapped. 

An advantage of simulations over experiments is that 
simulations have absolute knowledge of the system con- 
figuration. This allows tests of the scaling argument to 
be carried out, by testing its intermediate steps. It was 
found that Vdi ss /V ~ vCa and u ~ \[Q ~ -s/Ca. In order 
to obtain Vdiss and u it is necessary to determine a cut-off 
dtr for throat dissipations, with throats having d > d tI 
belonging to the dissipative volume. dt r is determined 
by increasing it to the point where the total dissipation 
contained above the threshold begins to decrease rapidly. 

The values for u and Vdiss obtained in this way are 
plotted for three values of Ca in Fig. [6j They scale as 
they are supposed to, in support of the scaling argument. 

In three dimensions, the perimeter of a stagnant cluster 
is a two-dimensional surface and the dissipative volume 
consists of sheets of flow surrounding these perimeter sur- 
faces. Postulates 1 and 2 do not need to be modified by 
this. Postulate 3 changes to the following: the sheets of 




Ca 



FIG. 6. (Color online) Scaling of mean velocity u (trian- 
gles) and dissipative volume Vdiss (circles) with Ca. Power 
laws (straight lines) are added with exponents 0.5. The units 
on the y-axis refers to the percentage of dissipative to total 
volume; velocities have arbitrary units. All simulations are 
within the stagnant cluster state, with S = 0.20 and M = 1. 



flow have a thickness which does not scale with Ca and 
the stagnant clusters obey the scaling I* ~ l x ~ l y ~ l z . 
Given these modifications the rest of the scaling argu- 
ment survives unaltered for three dimensions. The result 
P = 7 = 1/2 is therefore expected to be valid for both 
two and three-dimensional systems. 

An argument for a foam state exponent is not likely 
to share any similarities with the argument for the stag- 
nant cluster state. None of the three postulates apply to 
the foam state. There is no separation into two distinct 
volumes of low and high flow-rates, as can be seen from 
inspecting histograms of flow-rates from the two states. 
No regions appear to be stable with respect to either 
saturation distribution or flow-rates - flowing channels 
fluctuate across the entire system - which means that 
a balance criterion similar to postulate 2 is not applica- 
ble. The absence of stable structures makes it difficult to 
search for a characteristic length scale, while the power 
law distribution of flow-rates suggest that no character- 
istic length scale exists. 



VI. CONCLUSION 

Network simulations of two-phase flow in porous me- 
dia have been used to investigate the scaling of K e s with 
Ca under steady-state conditions with fixed saturation. 
Two values of M were used, corresponding to viscosity 
matched fluids and gas/liquid. Two states were identi- 
fied and named: the stagnant cluster state and the foam 
state. 
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The metastable stagnant cluster state scales with 7 = 
0.50 and 0.54 for viscosity matched fluids and gas/liquid 
respectively. An argument in support of 7 = 1/2 has 
been given. The argument applies to both the stagnant 
cluster state and the Oslo experiment. More work is 
needed to explore the scaling in the foam state. 

The Oslo experiment has been interpreted in terms of 
the distinction between forced and free fractional flow. 
Their low values of Ca and observed dependence between 
S and Ca is indication that they observe forced fractional 
flow, which indicates that their system is governed by the 
boundary conditions at the inlets. The Bozeman exper- 
iment, on the other hand, may have observed the foam 



state. 

Upper and lower bounds have been established for the 
power law. At Ca < Cai a spanning cluster transition 
prevents power law scaling for the case of free fractional 
flow, though the Oslo experiment indicates that scaling 
persists under conditions of forced fractional flow. At 
Ca > Ca up a viscous transition results in K e g ~ 1. The 
spanning cluster transition and viscous transition appear 
to be first and second order, respectively. 

Discussions with E. G. Flekk0y, A. Hansen, K. J. 
Mafoy, S. Sinha and K. T. Tallakstad are gratefully ac- 
knowledged. H. A. Knudscn is fondly remembered. 
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